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Superconducting phase transitions in strongly type-II superconductors in the Pauli paramagnetic 
limit are considered within the framework of the Gorkov-Ginzburg-Landau approach in the lowest 
Landau level approximation for both s and d-wave electron pairing. Simple analytical expressions 
for the quadratic and quartic coefficients in the order parameter expansion of the superconducting 
free energy are derived without relying on gradient or wavenumber expansions. The existence of 
a changeover from continuos to discontinuos superconducting phase transitions predicted to occur 
in the clean limit is shown to depend only on the dimensionality of the underlying electronic band 
structure. Such a changeover can take place in the quasi 2D regime below a critical value of a 3D-2D 
crossover parameter. 

PACS numbers: 74.20Rp, 74.25. Dw, 74.25.Bt, 74.78.-w 
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I. INTRODUCTION 

In recent years a number of unusual superconduct- 
ing (SC) states have been discovered in diflFerent new 
materialsi. Most of these materials are strongly type- 
II superconductors, possessing highly anisotropic or even 
quasi-two-dimensional (2D) electronic structures. Of 
special interest in the present paper are SC materials 
showing peculiar clean-limit features at high magnetic 
fields and low temperatures, notably the recently discov- 
ered family of heavy-fermion compounds CeRIns (R = 
Rh, Ir and Co)2>, and some of the organic charge transfer 
salts of the type (BEDT-TTF)2X3,4. 

The heavy-fermion compound CeCoIns, for example, 
which is believed to be an unconventional ( rf-wave) 
superconductor"" similar to the high-Tc cuprates, exhibits 
the highest Tc (~ 2.3 K) among the Ce-based heavy- 
Fermion compounds. This material is characterized by 
exceptionally strong Pauli paramagnetic pair-breaking^jii 
due to its extremely large electron effective mass and 
small Fermi velocity, which could lead to discontinu- 
ous (first-order) SC phase transitions at sufficiently high 
magnetic fields^f^fi^. 

Recently Bianchi et.alr^ have observed a dramatic 
changeover of the second-order SC phase transition to a 
first-order transition in specific heat measurements per- 
formed on this material as the magnetic field is increased 
above some critical values for both parallel and perpen- 
dicular field orientations with respect to the easy con- 
ducting planes. Similar effect has been very recently 
observed by Lortz et.al}^ in the nearly 2D organic su- 
perconductor K-(BEDT-TTF)2Cu(NCS)2, but only for 
magnetic field orientation parallel to the superconducting 
layers, where the orbital (diamagnetic) pair-breaking is 
completely suppressed. Under these conditions the usual 
(uniform) SC state is expected to be unstable with re- 
spect to formation of a nonuniform SC state, predicted 
more than 40 years ago by Fulde and Ferreli^,, and by 
Larkin and Ovchinikov^^ (FFLO). The corresponding SC 
order parameter is spatially-modulated along the field 



direction with a characteristic wavenumber, q , whose 
kinetic energy cost is compensated by the Pauli pair- 
breaking energy. The critical temperature, Tjfio, for the 
appearance of the FFLO phase is found to equal 0.56Tc. 
At the corresponding tricritical point the normal, the uni- 
form and non- uniform SC phases are all met. 

The possibility of a changeover to first-order transi- 
tions can be effectively investigated within the Ginzburg- 
Landau (GL) theory of superconductivity since for the 
uniform SC phase (i.e. for q — 0) the coefficient (usu- 
ally denoted by /3) of the quartic term in the GL expan- 
sion changes sign at a temperature T* , which coincides 
with Tffio^^. The identity of T* with Tffio is peculiar to 
the clean limit of a superconductor with no orbital pair- 
breaking. In conventional s-wave superconductors elec- 
tron scattering by non-magnetic impurities shifts T///o 
below the critical temperature r*i^j allowing discontin- 
uous phase transitions at temperatures Tffio <T <T* , 
since (following Anderson's theorem) j3 is not influenced 
by nonmagnetic impurities. In superconductors with un- 
conventional electron pairing, where (3 is strongly influ- 
enced by non-magnetic impurity scattering, the situation 
is reversed, i.e. T* < Tffio- 

The interplay between orbital and spin depairing in a 
pure s-wave isotropic 3D superconductor was flrst dis- 
cussed by Gruenberg and Gunther— , who conjectured 
(i.e. without presenting any result for the coefficient P 
) that for T < 0.56Tc the N-SC transition is of the sec- 
ond order whereas at lower field there should be a first- 
order transition to a uniform SC phase. Houzet and 
BuzdirtiSi have essentially confirmed this picture by ex- 
ploiting order-parameter and gradient expansions in the 
GL theory to find that T* < Tffio so that at tempera- 
tures T* < T < Tffio, there are second order transitions 
to either the LO or FF phase. It should be noted, how- 
ever, that the orbital effect was treated there by using 
gradient expansions, which is a valid approximation only 
at very low magnetic fields. 

In contrast to all the works outlined above, Adachi and 
Ikeda have recently foundi§' that, in a clean, d-wave, 2D 
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(layered) superconductor, the orbital effect always shifts 
Tffio below T* . In this work the authors have used or- 
der parameter expansion in the Gorkov Green's function 
approach to the GL theory up to six order, avoiding the 
restrictions of gradient expansion by exploiting the low- 
est Landau level (LLL) approximation for the conden- 
sate of Cooper-pairs. Accounting for impurity scattering 
destroys the FFLO phase and, in contrast to the pure 
paramagnetic situation, somewhat reduces T*. The ef- 
fect of SC thermal fluctuations was found in this work 
to broaden the discontinuous mean-field transition at T* 
into a crossover. The reliance on (FFLO) wavenumber 
expansion and on extensive numerical computations in 
this work has saved formidable analytical efforts, leaving 
however, interesting questions unanswered. In particu- 
lar, the origin of the relative shift of Tf fio below T* by 
the orbital effect found in this work, in contrast to all the 
other works, remains unknown. 



In the present paper we develop a formalism based 
on order parameter expansion within the Gorkov the- 
ory for a strongly type-II superconductor, with both s- 
and d^2_y2--wave electron pairing at high magnetic fields, 
which is sufficiently simple to yield useful analytical ex- 
pressions for the SC free energy to any desired order in 
the expansion. The fundamental interplay between spin 
induced paramagnetic and orbital diamagnetic effects at 
an arbitrary magnetic field is studied, within a model of 
anisotropic electron systems covering the entire 3D-2D 
crossover range, without relying on gradient or wavenum- 
ber expansions. These advantages enable us to shed new 
light on the yet undecided debate concerning the order 
of the SC phase transitions in the presence of strong Zee- 
man spin splitting, and to push our investigation into 
the unexplored region of very low temperatures, where 
quantum magnetic oscillations have been shown to be 
observable in the heavy fermion compounds2£i2i. 



Specifically, it is found that the relevant parameter 
controlling the relative position of Tf fio with respect to 
T* is the dimensionality of the electronic orbital motion 
in the crystal lattice, through its influence on the orbital 
(diamagnetic) pair-breaking effect. For a 3D Fermi sur- 
face (isotropic or anisotropic), where the electron motion 
along the magnetic field direction reduces the cyclotron 
kinetic energy, the shift of T* to low temperatures is 
larger than that of Tf fio- In this case the kinetic energy of 
Cooper-pairs associated with their motion along the field 
can compensate the spin-splitting effect, and thus leading 
to an increase of f3 and disappearance of the first-order 
transition. The corresponding phase diagram is similar 
to that suggested in Ref.^^, where the N-SC transition 
is of second order, whereas the transition between non- 
uniform and uniform (along the field) SC states is of first 
order. In the quasi-2D limit (i.e. for quasi-cylindrical 
Fermi surfaces) the enhanced orbital pair-breaking shifts 
Tffio below T* , in agreement with Adachi and Ikedcfci^. 



II. ORDER PARAMETER EXPANSION IN THE 
PRESENCE OF SPIN-SPLITT LANDAU LEVELS 



Our starting point is an effective BCS-like Hamiltonian 
with a c?2,2_j,2-wave pairing interaction similar to that ex- 
ploited, e.g. by Agterberg and Yang^S. The conventional 
s-wave situation can be similarly worked out and so will 
not be presented in detail here. The thermodynamical 
potential (per unit volume) for the corresponding d-wave 
superconductor, as expanded in the order parameter with 
nonlocal normal electron kernels, may be written as: 



m=l 



m 



-r!2™{A(R, r)} 



(1) 



where fl2m, {A (R, r)} is a functional of the SC order 
parameter, A(R, r), having a power-low dependence 
~ |Aop™ on the global amplitude, Ao , of the order 
parameter, and is a BCS coupling constant (given in 
units of energy X volume). The corresponding d-wave or- 
der parameter depends on both the center of mass (R) 
and relative (r) coordinates of a condensate of electron 
pairs: A(R, r) ~ A(R)(^(r). It should be determined 
self-consistently from the corresponding pair-correlation 
functions. Only stationary solutions are considered, ne- 
glecting quantum and thermal fluctuations. In addition 
the order parameter in the mean fleld approximation is 
selected as a hexagonal vortex lattice. Actually this as- 
sumption is not very important since the second order 
term in the order parameter expansion does not depend 
on the vortex lattice structure whereas the lattice struc- 
ture dependence of the quartic term is very weak (for 
reviews see^^^). 

For the underlying system of normal electrons we 
assume a simple model of quadratic energy disper- 
sion e (kx,ky,kz) = Y? (fc^ -I- /2m* + h^kll2ml and 
anisotropic effective mass tensor: m* < m*. A quasi- 
2D situation is characterized by a sufficiently large 
anisotropy parameter Xa = \/m*/m*, corresponding to 
an elongated Fermi surface with a Fermi momentum 
kp and Fermi energy ep = h'^kp/2m* , which is trun- 
cated by the Brillouin zone (BZ) face at kz^max = ir/d 
, where d is the lattice constant perpendicular to the 
easy planes. A parameter determining the dimensional- 
ity of the Fermi surface may be defined by: vo — 



, where £2, max = ^^k1 ^^^/2mZ, is the maximal value of 
the electron energy along the field. Thus, in the 2D limit, 

max <IC fci? , we have — > , while the system may be 
regarded 3D (isotropic or anisotropic ) if fcz,max — kp 
for which the Fermi surface is contained entirely within 
the first BZ, namely for = 1. 

At any order of the expansion, Eq.(IT]), the nonlocal 
electronic kernel of the corresponding functional (see 
e.g. Eqs. dS]),©, and Eqs. ([251) , (EH) ) consists of 
a product of to = 1,2,... pairs of normal electron 
Green's functions in a constant magnetic field, H — 
Hz (i.e. perpendicular to the easy conducting layers). 
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which are written in the form: G|| (Ri,R2,Wy) = , and the gauge invariant part can be calculated by the 
Gq^I (R2 — Ri, Wy) g (Ri, R2), where the gauge factor is well known expression^i 

given by 5(Ri,R2) = e^^^^'"'^^' "', an = \l ch/ eH 



ZTTajf J 2tt ^ fjL — Enfe^Ti + inwu + isign [lo,j) hi 

Here uji, = ttUbT {2v + 1) /h with v = Q, ±1, ... is Matzubara frequency, £n/c^| = hwc (n + 1/2 + a:^ — g/2) , Snk^i — 
u)c{n + 1/2 + + g/2), the spin-split normal electron energy levels, lOc = eH/m*-th.e in-plane electronic cyclotron 
frequency, x'^ — S/^k^ = 2m'uj ^ dimensionless longitudinal (parallel to the magnetic field) kinetic energy, bjf,g = 
eH/mo, is the Zeeman spin splitting energy, T- the impurity scattering relaxation rate, and /i = fiiUc {np + 1/2) w 
£f = h'^kp/2m,* is the chemical potential. The spatial variables are dimensionless in-plane (perpendicular to the 
magnetic field) coordinates, p — ^^^-R-u, ^ ^^^^ longitudinal coordinates: Zi = Hi ■ z , and Z2 = R2 • z. 



A. The quadratic term 

In the expansion, Eq.([T]), the second order term, which describes the SC condensation energy of spin-singlet electron- 
pairs, propagating from initial (i = 1) to final {i = 2) coordinates Ri ± is given by: 



^2 = ^-^1 d'Rid^R2t2 (Ri, R2) K2 (Ri, R2) = ^ - ^oA^ 

V Vo J V 



(3) 



where Vo — SLz is the volume of the system. The vertex part, r2, is a product of two order parameters multiplied 
by the gauge factors, g (R2, Ri), which are functions of the center of mass coordinates only, due to cacellation by the 
corresponding phase factors of the order parameters, namely: 

f 2 (Ri, R2) = 5* (Ri, R2) (R2, Ri) A (Ri) A* (R2) (4) 

The kernel K2 is a product of two translational invariant Green's functions, convoluted with the corresponding factors 
of the order parameters, which depend only on the relative pair coordinates, namely: 

^2(Rl,R2) f d''vi£v2^{vi)^*{v2) (5) 

X ^R2 — Ri + -^—^ — - , Lu^^ Gox ^Ri — R2 + — - , 

The factor of the order parameter which depends on the pair center of mass coordinates is written as2^: 

A(R) = c(Z)^e^^"'/20„(R^); 0„(R^) = e^^^-(^-5fr; a, = (6) 

n 

where c{Z) — coe*?^ is the Fulde-Ferrell modulation factor. Exploiting the fact that the kernel K2 (Ri,R2) depends 
only on the difference Ri — R2 , one may carry out the integration in Eq. ([3]) first over the in-plane mean coordinates 
Rj^= (R_L.i + R_L.2) /2 to get the following average vertex part2^ 

f 2) = ^ / f 2 (Rl, R2) d^R± = |co|2^e-'''/2-«<?(22-Zi) ^ ^2^-p-/2-^,{Z,-Z,) (7) 

' V J \'2n 

where Aq = ^ J d^R\A (R) p (Vo = SLz), and then integrate over the rest of the coordinates p = (R±,2 — R-^.i) /(^h 
and pz = (^2 - Zi) /an- 

Since, among other things, we are interested in the effect of quantum magnetic oscillations, we apply a tech- 
nique of exact summation over LLs suggested in Ref.^^. It is similar to the Poisson summation formula, which 
transforms the summation over LLs into summation over harmonics of the inverse magnetic field, and allows to 
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deal seperately with the uniform (quasi-classical) contribution and the various quantum corrections. This tech- 
nique can be briefly described as follows. Let us consider the integral representation of the Green's functions, 
\np — n — ± iuo\ ^ — dre^'^t"^"""^ **'^] and perform the summation over LLs using the well known identity, 
z"L„ (t) = (1 — z)^^exp (^jzjj with z = e*'"^ and t = Taking advantage of these relations the gauge 

invariant part of the Green's function for lo^ > can be transformed to: 

Gon (R. - = / ^e-^(-^-^) / dr ' , exp ( ) (8) 



2TTajjfuVc J 2tt Jq (1 - e-'^) V 4 1 

where lJi, = uj,y/ujc, T = T/ujc- For Ui, < one should replace r with — t (or Ui, with — |a;i,|). 

The scattering of electrons by non-magnetic impurities is taken into account here as a self-energy correction to 
the single electron Green's functions using the standard relaxation time approximation. Vertex corrections to the 
quadratic kernel if 2 (R-i,R2) (as well as to higher order ones), which are known to exactly cancel the self-energy 
insertions in the very weak magnetic field regime of convensional s-wave superconductors (see e.g.-^), are not so 
crucial in the strong magnetic field regime of both the s and d-wave situations investigated here, and will be therefore 
neglected in our calculations, as done, e.g. in Refsj ^^i^^ . In any event, for the high magnetic field and relatively 
clean superconductors considered here, the length scale, an, corresponding to the diamagnetic pair-breaking is much 
smaller than the electron mean free path vp/^, and the effect of impurity scattering is marginal. 

Utilizing this approximation we rewrite the kernel in the following form: 

K.{P,P.) ^ --^^E / d^^dz, [ / ^e--(--/^) (9) 

Jo Jo 



where a = (r_L,2 - rj_,i) /qh , ctz = (z2 - zi) /an, and 



J(Ti,T2,p) = y"dV,id'ri.2r (r2)(l-e-'"^) ' exp |^ 



{p + a/2f 1 + e"^ 
4 1 - e'^i 



In Eq.® we use the representation Lp{r) = 5 [z) f (rj_) , where / (rj^) = / /ke'*-'^''"^-* describes the two types 
of the electron pairing, the symmetric s-wave pairing with 

fsk = ^ (cos (kxd) + cos {kyd)) (11) 



and (i^2_y2-wave pairing. 



fdk = ^ (cos (k^d) - cos {kyd)) (12) 



producing nodes in the order parameter along the k^ = ±ky directions. The (5-dependence on z enables us to readily 
perfom the first two integrations in Eq. ([5]). 

Using the resulting expression for K2 {p, Pz), and the vertex function 2^ one can calculate the nontrivial coefficient 

Aq in Eq.Q by performing the integrals over p,and Pz 

Ao^ J d^pdp,K2 {p, Pz) e-^'/^-^^P. (13) 

with the other integrals incorporated in the kernel as appearing in Eq. ([9]) . It is convenient to perform the integration 
over p first since both the function J (n, r2,p)and the vertex part have a gaussian dependence on p which can be 
readily carried out with the result: 



J(ri,r2) = J d^pe-f"/^J{Ti,T2.^ 



^ 27rJ,(n,r2) 
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where 

Jp{ri,T2)^ J d^r^^,d^r^^2f {ri) f* {r2)e-^^" (15) 

and 7t = 2-e-'^i-e'^2 • should be noted here that the type of pairing influences the SC condensation energy 
through the functional dependence of Jp on the pairing function /k. 

Performing the straightforward calculation of Jp (ti, T2) for both functions one obtains, 

A,(n.T2) = i(l+e-J-"''")'==l (16) 



i.,(n,r,).i(l-.-J--)^i(2^) (17, 



where the last approximate step is obtained in the limit ^(P <C 1. This can be justified by noting that the scale of 
the function jj- is of the order unity whereas d (in units of the magnetic length) is much smaller than one. In the 
opposite limit: Jsp (n, T2) = Jdp (n, T2) = j. 

Thus, noting that the integration of ^0 over the center of mass coordinates yields just the total volume of the 
system, and performing the integration over the relative coordinate pz, 

^ 7 27r 7 27r 

27r 

one obtaines for a d-wave superconductor: 

QksT ■r-^ [ dkz 



A[^) ^ ""^^ , y / ^ / dr,dr2 (19) 

giTi[nF-^^(fe,+9/2)^+3+S„+ir]g-ir2[n^.-4^(fc,-g/2)^-g-ii2„-ir] 
X o . 

(2 - e^'^i - e*^2)'^ 

A similar expression can be derived for an s-wave superconductor. Below we will present only the final result for 
this case ( see Eqs. I^^ll^ . 

Eq. (jl9p is an exact representation for the coefficient of the quadratic term in the order parameter expansion, 
Eq|31 which includes low temperature quantum corrections and quantum magnetic oscillations. It can be written as 
a sum of contributions from poles at the 2D lattice: ti = 27rni and T2 — 27rn2 , with ni^2 = 0, 1, .... The dominant 
(zero harmonic) quasiclassical contribution arises from the pole at ni = n2 = whereas the quantum corrections are 
associated with the poles at ni = 712 ^ 0. It is easy to see that the oscillating terms correspond to the off-diagonal 
poles, ni 7^ 712, in the (ri, r2)-plane. 

In the present paper we are interested mainly in the quasiclassical contribution for which a further simplification can 
be achieved. Changing to new variables: T2 = po + f; ti = Po— §, and exploiting the expansion 2— e~*'^i — e*'^^ ~ Po"*''' 
near the "quasiclassical" pole ti = T2 = , one carries out the integral over t in Eq. (fT9|) to have: 

a]j , J Jo 

X 'np - e {kl + (g/2)')] ' e-P>--^"^''"+^'>'^^")] (20) 

Note that the lowest order expansion of the denominator in Eq. (fT9|) about ti = T2 = is kept under the entire 
range of integration since the important integration inerval is of the order r ^ pg ^ po ~ ^' ^'^^^ ^^^^^ that 

throughout this paper we assume that np ^ 1. 
It is convenient to rescale variables as 

u = ^/tlpPq; Xo = ^q; v= 4=^ (21) 
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and neglect the energy of an electron pair along the z-axis, (C?)^; with respect to Fermi energy, np. Perfoming the 
explicit summation over Matsubara frequencies one obtains, in terms of the new variables, the following result: 

Aid) 27:kBT r°° 1 - e-;;fe" / 2g \ ^(d) , . 

A-q' — N (0) \d—j=— — / du r-e cos\— —u\k3 2 \u,xo) (22) 



where = 3 (^^)^, ^^(0) is the electron density of states per spin at the Fermi energyC A^(0) = 2Trh^ " ~ )' 

ei^'^{u,xo)^ du{l~v^ycos{2xouv)e-^'(^-^'). (23) 
Jo 

A similar result is obtained for an s-wave superconductor. In this case ^ A^. = 1 and 

e^''(M,xo)= / dz^cos(2a;ouv)e-"'(i-'''). (24) 
Jo 

For 5 = xq = r = Eqs. ([22|) ([24 | reduces to the quadratic term derived by Helfand-Werthammer— . 

B. The quartic term 

The quartic term in the perturbation expansion, Eq.(IT]), which corresponds to a closed loop diagram with four 
vertices, is given by: 

r!^'^' = ^ / rf^Rld''R2rf^R3rf^R4r4(Rl,R2,R3,R4)i^4(Rl,R2,R3,R4), (25) 
1^0 J 

where the kernel, containing the gauge invariant factors of the four electron Green's functions, is: 
i'4(Ri,R2,R3,R4) - fcsT^ J dhid'^r2dh3dh4^{ri)ip* (r2)(^(r3)^* (r4) 



/ r2 — ri \ / r3 — r-? 

xGqi I R2-R1 H — .,uj„j Goi I R3 -R2 - ' ^ " 

X - R3 + , u;,^ Goi (ri - R4 - , j (26) 

and the vertex part: 

r4 (Ri, R2, R3, R4) = 9* (Ri, R2) 9 (R2, R3) 9* (R3, R4) 9 (R4, Ri) A (Ri) A* (R2) A (R3) A* (R4) (27) 

which consists of the gauge factors .g(Ri, Rj) and the order parameter values at the four center of mass positions for 
two electron pairs. 

Since the dependence of the order parameter on the relative pair coordinates is separable from that of the center 
of mass coordinates, the latter dependence is selected to have the usual Abrikosov lattice structure, 

A(R) = Co e''''^e-Kl"l'-«') ^ e*9,.«-9S/4^ (28) 

n=0,±l,±2,... 

with g„ = 21^71/ ax and u = X + iY, . To simplfy the calculation of the vertex part we exploit several assumptions. 
Substituting Eq. ([28|) to Eq. ([27|) . one may keep only diagonal terms with g„i — g„2 = <ln3 = <lnA = P , since all off- 
diagonal terms are small by the gaussian factor ^ exp — (g„4 — (?rii)^ — {<lnA — 9ni)^ ■ Furthermore, we may replace 
summation over p with an appropriate integration. Both of these assumptions are equivalent to neglecting particular 
vortex lattice structures, corresponding to replacement of the Abrikosov structure parameter, Pa, with which 
yields only a small error. 
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With the above assumptions the vertex part reduces to: 

r4(Ri,R2,R3,R4) = ^^^|co|V?(^^-^=+^-'-^^)e-^^'-il'"l%3[(«i-«3)^+K-<)^] (29) 

ZTT 

where pi = uz+i — ui. Since the dominant contribution to the quartic term arises from small propagation distances, 
\ui\ < 1— i^, one may expand the last exponential on the RHS of Eq. ip^)) . up to leading order, under the integrals 
over angular variables in Eq. (|25|) . Additional angular dependence is due to the kernel, K4, through its dependence on 
the absolute values of linear combinations of "external" , R;+i — R;, and "internal" , rj+i — r/ (/ = 1, .., 4), coordinates 
(see Eq. [26l) . Since the characteristic size of |ri+i — r;] '--^ d is much smaller than the scale of |Ri+i — Ri| '-^ qh , the 
dependence of the kernel on rj+i — r; (and consequently its dependence on the angular variables) may be neglected at 
large |R;+i — R;|. Therefore, the integration over angular variables in this region involves only the last exponential 
in Eq. (pg|) . resulting in: 

(eM("i-"3)^+W-<)^]^ « 1 + i((7ii-W3)' + (7/; -<)') = 1, 

since (uf) = {uiuD = , {I ^ k) 

whereas for small values of |Ri+i — R;| this exponential is always close to 1 and the remaining integration over angular 
variables can be perfomed in closed form (see below). Thus, one can approximate the vertex part by the following 
simple expression: 

r4 (Ri, R2, R3, R4) = ^^|co|^e^«(^^-^^+^^-^^)e-3 2: Ip,P (30) 

ZTT 

which depends only on nearest neighboring coordinates. 

Making use of Eq. pop . the remaining calculation of the quartic term is similar to that used for the quadratic term, 
but considerably massier. Below we present only an outline of the derivation. Since integrations over Zi are trivial we 
shall use from now on only 2D vector notations with integrations over Zi written explicitly. 

Combining Eqs. (|25|) .([26 |) . ([30l) our starting expression for the quartic term is given by: 

Js,d) _ ksT f _L X"* ax\^ 



X 



n^^^g-ifc,,i(Z2-Zi)glfc,,2(23-22)g-ifc,,3(Z4-Z3)g»fc.,4(2l-24)gi9(Zl-22+23-24) 
2n 

i=l 



X J rf2ridWr3d'r4/(ri)r (r2)/(r3)r (r4) X e4(ri,r2,r3,r4;{fc,,J;cc.,) (31) 
where the function O4 (ri, r2, r3, r^; {k^^i} ; uj^) includes integration over all electron pair coordinates: 
e4(ri,r2,r3,r4;{fc^,J;Lj,) = / d2p_^^2j^^^2j^^^2j^^g-i E IftP x 



drie 



exp 



V 4 1 



-X2+g~{-iujj^+iT 



exp 



-ff-23 l+e'"2 
4 l~e"2 



__expf^i±^) poo _ exp f 

Jo 1 - e-'^4 ^ ^ 

Here the coordinates, R^.i+i, in Ea. p2p are the linear combinations of pi = R(+i — R; and rji — r;+i — r; : 

R12 = Pi + -jVi\ R23 = P2 - 2^72 

R34 = P3 + ^?73; R41 = PA - ^V4 (33) 
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The gaussian integration over pi reduces Eq. ([32)1 to: 

Jo ^ 



ax 



(34) 



e . It should be noted here that Eq.(f34|) has been obtained by exploiting the 



where x — A: — e"^ — e~"^ — e*'^^ 
fact that the dominant contributions to the integrals originate in the regions where <C 1. 

Furthermoe, noting that in the above equation the rji- and fc^-dependences are factorized, one can perfom the 
integrations over both sets of variables separetely. For a d-wave superconductor we obtain: 



j dh,d\2dh3d\if (n) /* (ra) / (ra) /* (r4) x exp 



-F- ('71 -m + m- mf 

ox 



43 ( 



9 \ 4 



2 _ 9_(P_ 
~ 16^ 



(35) 



where the last appoximation is valid under the same conditions discussed in the derivation of the quadratic term. 
Thus the quartic term is transformed to: 



ksT ax^/ir . ,4 
[huoc) (27r) 



16 



dTldT2dTy,dT/^ 



(4-e 



^-i(ri-T2+T3-r4)(nF-4^fc^-^^(f )^)^_(ri+T2+r3+r4)(w„+r)gi(Ti+r2+r3+r4)(g+^^fc,g) ^^gg-j 

where an additional integration over C = ('''1 ^ 7'2 + T3 — T4) /2 for small Ti can be performed. Rescaling variables as 



Tl + T2+T3 + Ti 



2 ' 2 

and summing up over v one obtain the final result for the quartic term 



s = — :— (T3 - Ti) ; 



= B^^^ / dg 



1 — e 



sinhf2£^p 



e -y^^ cos {2 Qgo)ef'' {g,q) 



where S'^''^ = cf"^ Bg with 



Sd) _ _3_ fkpd] 
'4 16 V 51- J 



(e, 9) = dv (1 - v^)' e-i^'i'-^") cos (2^;g^,) (^^' dse-''(''^') 



(37) 



(38) 



The result for an s-wave superconductor can be obtained 
from Eq. ([57|) by replacing the factor (l — v"^) in the def- 
inition of qII^"^ and the factor 04^^^ in the normalization 
coefficient B^'^"> with unity. The s-wave quartic term for 
zero spin splitting is equivalent to that obtained in Ref.^^. 



III. RESULTS AND DISCUSSION 

The analysis presented in the previous sections enables 
us to write a GL-like expansion of the SC contribution 
to the thermodynamic potential for an s or dx2^y2- wave 
pairing up to second order in Aq: 



(39) 
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For the quadratic term we have: 



{t, b, q) 



1 



dp- 



1-e '(^o) 



smh ( 



dvd':^''^^ (") cos (qpv/xa) 
(1 - v') p^l2a\\ 



/o 
exp 



(40) 



where A = N{^)V , {Td) = hvp/TrkBTD , with To 
the Debye teraperature, and vp — ^y 2e p / m* -the inplane 
Fermi velocity, <; (T) = hvp /i^kBT- the thermal mean- 
free path, rp = \/1npaH- the electronic cyclotron radius 
at the Fermi energy, and I is the mean-free path due to 
impurity scattering. 

The differences between s-wave and d-wave SCs are 
given by 4'^ = 1, 4"^' = 3(^)^ and dS^^ [v) = 1, 

The quartic term has a similar structure: 



dp- 



_ 2p 

1-e ^(^d) 



e-^Pl' cos 



sinh(^) 



rp 



(41) 



e 



{s4) 



dv-d^^'"^^ (v) cos [qpv/xa) 



exp[- (l-.2)pV4ay ^ 



dse-^'i'-^') 



mjihBT (s) ^ ^ u) ^ 



)(<*) ^,^ - 



where nn 

^ (^)', (t;) = 1, and <^ («) = (1 - « ; • 

On the basis of the above formulas we discuss be- 
low the H-T phase diagram for different values of the 
relevant parameters. Three independent dimensionless 
parameters: {2aH /rp) g , 2a///<r(T) and qan/Xa, con- 
trol the basic integrals in these equations. The first two 
parameters measure the strength of the spin and ther- 
mal pair-breaking mechanisms, respectively, relative to 
the orbital (diamagnctic) depairing. The third parame- 
ter determines the relative strength of the compensating 
FFLO mechanism. The value of the spin pair-breaking 
parameter, a = g {^au /tp) fj^^jorh , where i?°20 ^^e 
upper critical field at T = in the absence of spin pair- 
breaking, is related to the well known Maki parameter^, 

"Af = ^''{[lekBT^r ' ^ l-lttM- Here T^q is the 
transition temperature at zero magnetic field. 

As we shall show below, the situation Tffio > T* , 
where T* is the temperature at which /? {t, bc2,q — 0) — 
, is realized in 3D systems (corresponding to t^o = 1), re- 
gardless of the spin-splitting strength and the type of 
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FIG. 1: Phase transition lines for a 3D system with s-wave 
pairing : N-SC transitions for order parameters with (solid 
line) or without (doted line) FFLO modulation, and FFLO- 
BCS transitions obtained for GL free energy with (dashed 
line) or without (doted-dashed line) quartic correction. The 
N-SC transition is of second order whereas the FFLO-BCS 
transition is of first order. The value of the spin splitting 
parameter is a — 1.8. 
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FIG. 2: The GL coefficients (in arbitrary units ) P {t,b,q = 0) 
(dashed lines), ''"'■'^^'j'^^-' (solid lines) ,^^^%|t- 
as functions of the parameter |^ for (a) vq 
vo = 1. 



— (doted lines) 
= .4 and (b) 



electron pairing. A typical phase diagram is shown in 
Fig. 1 for s-wave pairing and spin pair-breaking param- 
eter (7 = 3. 

As long as r > T* (so that /3{t,bc2,q^ 0) > 0) the 
normal to SC (N-SC) phase transition is of second-order 
and the (reduced) critical field, bc2 (t), can be determined 
as the maximal value of 6 = H/H°2q obtained from the 
equation a (t, b,q) = for all values of q , at the (reduced) 
temperature t = T/Tcq. The solution of this equation 

for q = yields a transition line, b = 6^2' (^)i ignoring 
the possibility of a FFLO state. The tricritical point, 
Tffio, is defined as the maximal temperature at which 

6c2 (t) - ^^'^ 



> (*) ■ It can alternatively be determined from 



the equation 



da{t,bc2,q—0) 



— , which is equivalent to the 



condition for vanishing of the coefficient of 
gradient expansion of the SC free energyi^.. 
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values the coefficient (3 (t,bc2, q) at optimal q is always 
positive and the N-SC transition is of the second order 
at arbitrarily low temperature. 

The transition within the SC region from the nonuni- 
form (FFLO) to uniform (BCS) phase at T > T* can 
not be obtained just by analyzing the quadratic term 
a it, b, q) since the SC order parameter is finite there. It 
can be obtained by minimizing the SC free energy (in- 
cluding both quadratic and quartic terms) with respect 
to the modulation wave number q. Neglecting the sixth 
and higher order terms in the expansion, the correspond- 
ing (standard) GL free energy, Q (q) ~ ~9 (a) jp, { (a) 
being the Heaviside step function), which has a single 
minimum aX q ^ for field near bc2 (see Fig. 3a), devel- 
opes a double- well structure (see Fig. 3b) as a function of 
q upon decreasing the field below 6c2 at a given tempera- 
ture T (due to the symmetry q <-> —q only positive values 
may be considered). One of these minima is always at 
q — 0, and it becomes energetically favorable at a critical 
field for a first order phase transition from the FFLO to 
the uniform BCS phase. The second (metastable) mini- 
mum at g ^ disappears completely upon further field 
decrease (see Fig. 3c). 



\ » -0.025 

\ / / 

A / 



FIG. 3: (a) The dependence of the GL coefficient a (orange) 
and the mean field SC free energy, —6 (a) (blue), on the 

modulation paramer, q= (|-) yj^^ {^o,h°^'^ ^ ^ 3D sys- 
tem (wo = 1) , a.t t = 0.4, b — 0.1142, i.e. near the tricritical 
point, just below the Normal-nonuniform SC transition. It 
is seen that a < in a small region around q — 0.6 where 
the SC free energy has a minimum. The value of the spin 
splitting parameter is cr = 3. (b): The same as in (a) but 
for a slightly lower field, b — 0.114, where a uniform (q = 0) 
metastable SC state is present, (c): The same as in (b) but 
for a slightly lower field, b = 0.1139, where a uniform {q = 0) 
equilibrium SC state is present, while a metastable SC state 
exists at g 7^ 0. 



For T < T* and sufficiently strong spin pair-breaking 
there can be a changeover to first-order SC transitions, 
but since /3 (t, 6c2, <? 7^ 0) > (3 {t,bc2,q — 0) (see Fig. 2 ), 
the segment of the 6c2 (i)-line with first order transitions 
arises only at very low temperatures. For moderate a 




FIG. 4: The GL coefficients, a (green) and /3 (blue), and 
the mean field SC free energy (orange) as functions of the 
modulation parameter, q, in a 3D system {vo = 1) at a rel- 
atively low temperature t — .25 and decreasing field values 
(a) b — .118 and (b) 6 = .117. Note the vanishing of /3 inside 
the region where q < , around which the used approxima- 
tion, —9 (q) 2§^y for f2(g) breaks down (dashed sector of the 
orange line ). 
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At temperatures T below T* the first two terms in 
the expansion of the thermodynamic potential are not 
sufficient to correctly describe the uniform SC state since 
for negative /3- values the scale of the SC free energy is 
determined by the sixth order term. In contrast, the free 
energy of the nonuniform state, where /3 (t, 6, 5 7^ 0) > 0, 
can be obtained from the stantard GL functional (with 
the assumption that the contribution of the sixth order 
term is small compared to that of the quartic term) . The 
characteristic g-dependences of the GL coefficients, a and 
/3, and the mean field free energy —9 {a) 2p^0 ; < 
are illustrated in Fig. 4. Whereas at high fields (Fig. 
4a) the minimum of the SC energy occurs in a region 
where > 0, at lower fields (see Fig. 4b) it approaches 
the expanding temperature domain of negative /3. Thus, 
even for moderate spin splitting and low temperature the 
transition line from the nonuniform to uniform SC state 
cannot be determined without knowing the sixth-order 
term. It is clear, however, that this transition is of the 
first order. 

It should be noted that if one attempts to deter- 
mine the FFLO-BCS phase boundary from the equation 
'^°'^*dq^~^^ = it will greatly overestimate the size of the 
FFLO phase as compared to that obtained by minimiz- 
ing — ^ (see Fig. 1). This remarkable difference is due 
to the strong g^-dependence of the quartic coefficient (3 
(see Fig. 2). 

The suppression of the orbital effect in the considered 
3D systems, with ellipsoidal Fermi surfaces contained en- 
tirely within the BZ, is due to the factor 1 — appearing 
in the Gaussian exponents of Eqs. (HO]) , (HH) . The re- 
covery of this effect in quasi-2D systems with truncated 
ellipsoidal Fermi surface, where Vq < 1, can reverse the 
relation between Tf/io and T* . Fig. 2b, where the GL 
coefficients are shown for a ~ 1.8, vq = .4 and s-wave 
pairing, illustrates the situation with Tffio < T*, which 
occurs for all values of vq below a critical dimensionality 
vo.cr ~ 0.44 (see Fig. 5), and depends only weakly on 
the spin-splitting parameter a. 

The corresponding phase diagram (see Fig. 6 ) for 
Vq below this crossing point is quite different from that 
found for the 3D systems shown in Fig.l. First of all, 
since /3 < , one may use Eqs. (1^0]) . ([1T|) to determine the 
phase diagram only under the assumption that the sixth 
order coefficient 7 is positive (see RefjiS). In this case a 
discontinuous SC transition occurs at n (Ag) = with 

Aq = (^^4^) and a = > , and the corresponding 

critical field, bc2{t) , should be larger than 6^2 (i) j ob- 
tained from the equation a(q = 0) 0. Thus, at a tem- 
perature below T* , the N-SC phase boundary includes a 
segment of first order transitions, which may end at zero 
temperature, or at a finite temperature, depending on 
the spin-splitting strength. This dependence appears be- 
cause of the competition between the decreasing explicit 
dependence of f3 on decreasing temperature and its in- 
creasing implicit dependence through q (T) at the FFLO 



H 



0.5 



FIG. 5; Solutions |^ cx T of the equation P {q = 0) 
responding to T* (dashed line), and the equation 
0, corresponding to T//io (solid line) vs the dimensionality 
crossover parameter vq for a = 2.5. 



0, cor- 

da(q^O) _ 




FIG. 6: Schematic phase diagram for a quasi 2D system. The 
shaded area corresponds to a nonuniform SC (FFLO) phase, 
the dashed line corresponds to a {t,b, q = 0) = , and the 
dotted-dashed Une can be obtained from -^^^a {t, b,q = 0) — 0. 



state. The boundary between the BCS and FFLO states 
should be determined by minimizing the free energy, Eq. 
([39)) . with respect to q. This may be restricted to the 
explicit dependence on q since the order parameter is de- 
termined by: — 0. Consequently the positive sign of 

'^^''*dq2~^^ (see Fig. 2b) results in partial cancellation of 

the leading contribution to , which is proportional to 

^^"^ negative in the FFLO part of the phase 
diagram. Moreover, since for the discontinuous transition 
the order parameter is finite just below the transition the 
higher order terms in Aq (Eq. ((39|) ) should be taken into 
account. As a result Tffio should be smaller than T^fio' 

the temperature obtained from the equation 
= , as schematically shown in Fig. 6. 
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IV. CONCLUSIONS 

It is shown that the expected changeover to first-order 
SC transitions in clean, strongly type-II superconductors 
in the Pauli paramagnetic limit can take place only in ma- 
terials with quasi-cylindrical Fermi surfaces, regardless of 
the type of the electron (s or d-wave) pairing interaction 
which leads to superconductivity. This finding clarifies 
the confusing current literature on this topioiifi^fi^. 

The observation of such a changeover in the heavy 
fermion compound CeCoIns for magnetic field orienta- 
tion perpendicular to the easy conducting plane^i is con- 
sistent with the quasi-2D character of its electronic band 



structure^^. The interesting situation of a 2D supercon- 
ductor under a magnetic field parallel to the conducting 
plane, for which a changeover to discontinuous SC tran- 
sitions was reported very recently^^, is more subtle since 
the vanishingly small cyclotron frequency characterizing 
this case does not allow utilization of the Landau orbitals 
approach employed here (for a recent review see, e.g.—). 
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